fish_traits_sum <- fish_traits %>%mutate(across(where(is.numeric), round, 2))## Display the table ----htmltools::tagList(DT::datatable(fish_traits_sum))
traits correlation
Code
M <-cor(numeric_traits[, c(-1)])ggcorrplot::ggcorrplot(M, hc.order =TRUE, type ="lower",lab =TRUE, tl.cex =9, lab_size =3)
Code
# list of species sp_names <-c(rownames(fish_traits))# taxonomic_familiestaxonomic_families <- sp_names %>%as.data.frame() %>%`colnames<-`("species") %>%mutate(family =case_when( species %in%c("Benthosema_glaciale","Ceratoscopelus_maderensis","Diaphus_metopoclampus","Lampanyctus_ater","Lampanyctus_crocodilus","Lampanyctus_macdonaldi","Lobianchia_gemellarii","Myctophum_punctatum","Notoscopelus_bolini","Notoscopelus_kroyeri","Bolinichthys_supralateralis" ) ~"Myctophidae", species %in%c("Borostomias_antarcticus","Chauliodus_sloani","Malacosteus_niger","Melanostomias_bartonbeani","Stomias_boa" ) ~"Stomiidae", species %in%c("Holtbyrnia_anomala","Holtbyrnia_macrops","Maulisia_argipalla","Maulisia_mauli","Maulisia_microlepis","Normichthys_operosus","Searsia_koefoedi","Sagamichthys_schnakenbecki" ) ~"Platytroctidae", species %in%c("Sigmops_bathyphilus","Gonostoma_elongatum") ~"Gonostomatidae", species %in%c("Argyropelecus_hemigymnus","Maurolicus_muelleri","Argyropelecus_olfersii" ) ~"Sternoptychidae", species =="Anoplogaster_cornuta"~"Anoplogastridae", species %in%c("Arctozenus_risso", "Paralepis_coregonoides") ~"Paralepididae", species =="Bathylagus_euryops"~"Bathylagidae", species =="Cyclothone"~"Gonostomatidae", species =="Derichthys_serpentinus"~"Derichthyidae", species =="Eurypharynx_pelecanoides"~"Eurypharyngidae", species =="Evermannella_balbo"~"Evermannellidae", species =="Lestidiops_sphyrenoides"~"Lestidiidae", species =="Melanostigma_atlanticum"~"Zoarcidae", species %in%c("Photostylus_pycnopterus","Xenodermichthys_copei") ~"Alepocephalidae", species =="Serrivomer_beanii"~"Serrivomeridae" ) )%>%mutate(order =case_when( family =="Myctophidae"~"Myctophiformes", family %in%c("Stomiidae","Gonostomatidae", "Sternoptychidae") ~"Stomiiformes", family %in%c("Platytroctidae","Alepocephalidae") ~"Alepocephaliformes", family =="Anoplogastridae"~"Trachichthyiformes", family %in%c("Paralepididae","Evermannellidae","Lestidiidae") ~"Aulopiformes", family =="Bathylagidae"~"Argentiniformes", family %in%c("Derichthyidae","Serrivomeridae") ~"Anguilliformes", family =="Eurypharyngidae"~"Saccopharyngiformes", family =="Zoarcidae"~"Perciformes", ) )
## Summary of the assemblages * species data.frame ----asb_sp_fish_summ <- mFD::asb.sp.summary(asb_sp_w = depth_fish_biomass)asb_sp_fish_occ <- asb_sp_fish_summ$"asb_sp_occ"htmltools::tagList(DT::datatable(asb_sp_fish_occ))
2.2 Computing distances between species based on functional traits
We have non-continuous traits so we use the Gower distance(metric = “gower”) as this method allows traits weighting.
2.3 Building functional spaces and chosing the best one
2.3.1 Computing several multimensional functional spaces and assessing their quality
mFD evaluates the quality of PCoA-based multidimensional spaces according to the deviation between trait-based distances and distances in the functional space (extension of Maire et al. (2015) framework).
2.3.3 Testing the correlation between functional axes and traits
Code
sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"# As we have 26 traits we have to split the df to see correlation between functional axes and traits # first set ----fish_traits_1 <- fish_traits%>%select(1:9)fish_tr_faxes <- mFD::traits.faxes.cor(sp_tr = fish_traits_1, sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], plot = T)## Print traits with significant effect ----fish_tr_faxes$"tr_faxes_stat"[which(fish_tr_faxes$"tr_faxes_stat"$"p.value"<0.05), ]
trait axis test stat value p.value
1 eye_size PC1 Linear Model r2 0.160 0.0095
2 eye_size PC2 Linear Model r2 0.191 0.0043
3 eye_size PC3 Linear Model r2 0.105 0.0388
4 eye_size PC4 Linear Model r2 0.160 0.0096
5 orbital_length PC1 Linear Model r2 0.414 0.0000
10 gill_outflow PC2 Linear Model r2 0.518 0.0000
13 oral_gape_surface PC1 Linear Model r2 0.126 0.0226
14 oral_gape_surface PC2 Linear Model r2 0.522 0.0000
18 oral_gape_shape PC2 Linear Model r2 0.160 0.0096
24 oral_gape_position PC4 Linear Model r2 0.245 0.0010
26 lower_jaw_length PC2 Linear Model r2 0.695 0.0000
29 head_length PC1 Linear Model r2 0.337 0.0001
30 head_length PC2 Linear Model r2 0.367 0.0000
34 body_depth PC2 Linear Model r2 0.353 0.0000
35 body_depth PC3 Linear Model r2 0.199 0.0034
Code
## Plot ----fish_tr_faxes$"tr_faxes_plot"
Code
# second set ----fish_traits_2 <- fish_traits%>%select(10:18)fish_tr_faxes_2 <- mFD::traits.faxes.cor(sp_tr = fish_traits_2, sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], plot = T)## Print traits with significant effect ----fish_tr_faxes_2$"tr_faxes_stat"[which(fish_tr_faxes_2$"tr_faxes_stat"$"p.value"<0.05), ]
trait axis test stat value p.value
2 pectoral_fin_position PC2 Linear Model r2 0.155 0.0108
5 pectoral_fin_insertion PC1 Linear Model r2 0.291 0.0003
6 pectoral_fin_insertion PC2 Linear Model r2 0.378 0.0000
9 transversal_shape PC1 Linear Model r2 0.122 0.0252
11 transversal_shape PC3 Linear Model r2 0.212 0.0024
14 caudal_throttle_width PC2 Linear Model r2 0.410 0.0000
19 dorsal_fin_insertion PC3 Linear Model r2 0.193 0.0041
23 eye_position PC3 Linear Model r2 0.201 0.0033
32 ventral_photophores PC4 Kruskal-Wallis eta2 0.590 0.0000
33 gland_head PC1 Kruskal-Wallis eta2 0.245 0.0011
Code
## Plot ----fish_tr_faxes_2$"tr_faxes_plot"
Code
# third set ----fish_traits_3 <- fish_traits%>%select(19:26)fish_tr_faxes_3 <- mFD::traits.faxes.cor(sp_tr = fish_traits_3, sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], plot = T)## Print traits with significant effect ----fish_tr_faxes_3$"tr_faxes_stat"[which(fish_tr_faxes_3$"tr_faxes_stat"$"p.value"<0.05), ]
Standard Effect Size (SES): to eliminate the influence of species richness on the functional diversity indices (Mouchet et al., 2010). Measures the deviation from the random expectation in standard deviation units
null model frequency: Randomize community data matrix abundances (here biomasss) within species (maintains species occurrence frequency)
FRic Functional Richness: the proportion of functional space filled by species of the studied assemblage, i.e. the volume inside the convex-hull shaping species. To compute FRic the number of species must be at least higher than the number of functional axis + 1.
FDis Functional Dispersion: the biomass weighted deviation of species traits values from the center of the functional space filled by the assemblage i.e. the biomass-weighted mean distance to the biomass-weighted mean trait values of the assemblage.
FDiv Functional Divergence: the proportion of the biomass supported by the species with the most extreme functional traits i.e. the ones located close to the edge of the convex-hull filled by the assemblage.
FEve Functional Evenness: the regularity of biomass distribution in the functional space using the Minimum Spanning Tree linking all species present in the assemblage.
calcul des indices avec package FD
correction distance matrix (obtenue avec Gower) = “none” or “lingoes” transformation avec “sqrt” ne permet pas de transformer la matrice en “euclidean”
randomisation des traits des espèces par couche de profondeur
Code
# Model parameters ----# Number of simulationsn_simulations <-999# correction method to use when the distance matrix cannot be represented in a Euclidean spacecorr_method <-"lingoes"# Depth layersdepth_layers <-rownames(depth_fish_biomass)# List of indices to be calculatedindices <-c("FRic", "FDis", "FDiv", "FEve")# Matrices for storing observed and simulated resultsdbFD_result_obs <-matrix(NA, nrow =length(depth_layers), ncol =length(indices),dimnames =list(depth_layers, indices))dbFD_result_sim <-array(NA, dim =c(length(depth_layers), n_simulations, length(indices)),dimnames =list(depth_layers, paste0("Sim.", 1:n_simulations), indices))dbFD_result_sim_list <-list()randomize_traits <-function(traits_mat) { randomized_mat <- traits_mat# Randomisation colonne par colonnefor (trait incolnames(traits_mat)) { trait_values <- traits_mat[[trait]]# Randomisation avec remise randomized_mat[[trait]] <-sample(trait_values, size =length(trait_values), replace =TRUE) }# Supprimer les niveaux vides des facteurs randomized_mat <-droplevels(randomized_mat)return(randomized_mat)}# Loop on each depth layer ----for (layer in depth_layers) {# Select the species present in the selected depth layer species_in_layer <-colnames(depth_fish_biomass)[depth_fish_biomass[layer, ] >0] traits_layer <- fish_traits[species_in_layer, , drop =FALSE] %>%select(where(~length(unique(.)) >1)) %>%droplevels()# If no variable lines remain, we move on to the next layerif (ncol(traits_layer) ==0) next# Correct formatting of biomass_layer with depth layer in rownames# Converts biomass data to matrix, which is required by FD::dbFD() biomass_layer <-matrix(as.numeric(depth_fish_biomass[layer, species_in_layer, drop =FALSE]),nrow =1, dimnames =list(layer, species_in_layer))# Simulations for each layerfor (sim in1:n_simulations) { randomized_traits <-randomize_traits(traits_layer) dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun =TRUE, m =4, a = biomass_layer, messages=F,corr = corr_method)# Storage of simulated values dbFD_result_sim[layer, sim, "FRic"] <- dbFD_sim_result$FRic dbFD_result_sim[layer, sim, "FDis"] <- dbFD_sim_result$FDis dbFD_result_sim[layer, sim, "FEve"] <- dbFD_sim_result$FEve dbFD_result_sim[layer, sim, "FDiv"] <- dbFD_sim_result$FDiv# Save simulated results for normality tests dbFD_result_sim_list <- dbFD_result_sim }}# ---- FIN ----
3.2. Statistic tests
Code
## Calcul observed values observed_values <- FD::dbFD(x = fish_traits, w.abun =TRUE, m =4, a = depth_fish_biomass,messages =FALSE, corr = corr_method) observed_values <- observed_values[c("FRic", "FDiv", "FDis", "FEve")] %>%as.data.frame() %>% tibble::rownames_to_column(var ="depth_layer") %>% tidyr::pivot_longer(!depth_layer, values_to ="observed_values", names_to ="indice")# Test normality ----## Arrange simulated values simulated_values <- dbFD_result_sim_list[, , ] %>%as.data.frame() %>% tibble::rownames_to_column(var ="depth_layer") %>% tidyr::pivot_longer(!depth_layer, values_to ="values", names_to ="indice") %>%mutate(indice =gsub("Sim\\.\\d+\\.", "", indice))# Statistic testsstat_test <- simulated_values %>%group_by(indice, depth_layer) %>%mutate(skewness_values = e1071::skewness(values),shapiro_p_values = rstatix::shapiro_test(values)$p.value ) %>%select(-values) %>%distinct() %>%# which simulated diversity indices have normal distributionmutate(normality=case_when(shapiro_p_values>0.05~"normal", shapiro_p_values<=0.05~"non_normal"))test_values_trans <- simulated_values %>%inner_join(stat_test %>%select(normality)) %>%filter(normality =="non_normal") %>%mutate(values_log =log10(values)) %>%group_by(indice, depth_layer) %>%mutate(skewness_values = e1071::skewness(values_log),shapiro_p_values = rstatix::shapiro_test(values_log)$p.value ) %>%select(depth_layer, indice, shapiro_p_values) %>%distinct() %>%# which simulated diversity indices have normal distributionmutate(normality =case_when( shapiro_p_values >0.05~"normal", shapiro_p_values <=0.05~"non_normal" ))values_percentile <- simulated_values %>%inner_join(stat_test %>%select(normality)) %>%filter(normality =="non_normal") %>%mutate(values_log =log10(values)) %>%group_by(depth_layer, indice) %>%summarise(quantile = scales::percent(c(0.025, 0.975)),percentile =quantile(values_log, c(0.025, 0.975))) %>%distinct() %>%inner_join(observed_values) %>%mutate(observed_values_trans=(log10(observed_values))) %>%select(-observed_values) %>%filter(indice=="FRic")# SES calcul ----## calcul mean and sd of simlauted values not transformsum_simulated_values <- simulated_values %>%# simulated values of indices normaly distributed (or not but for which the log doesn't change it)filter(indice%in%c("FDiv", "FDis", "FEve")) %>%group_by(indice, depth_layer) %>%mutate(meanNullFD =mean(values), sdNullFD =sd(values)) %>%select(depth_layer, indice, meanNullFD, sdNullFD) %>%distinct()## SES calcul not transformSES_calcul <- observed_values %>%inner_join(sum_simulated_values) %>%group_by(indice, depth_layer) %>%mutate(SES= (observed_values-meanNullFD)/sdNullFD) %>%filter(indice%in%c("FDis", "FEve", "FDiv"))## calcul mean and sd of simulated transformed values sum_simulated_values_log <- simulated_values %>%# simulated values of indices not normality distributed filter(indice=="FRic") %>%mutate(values_log =log10(values)) %>%group_by(depth_layer) %>%mutate(meanNullFD =mean(values_log), sdNullFD =sd(values_log)) %>%select(depth_layer, indice, meanNullFD, sdNullFD) %>%distinct()SES_calcul_log <- observed_values %>%inner_join(sum_simulated_values_log) %>%group_by(depth_layer) %>%mutate(SES= (log10(observed_values)-meanNullFD)/sdNullFD)SES_combined <- SES_calcul_log %>%full_join(SES_calcul) %>%select(-c(meanNullFD, sdNullFD))# plot ----SES_combined$depth_layer <-factor( SES_combined$depth_layer,levels =c("Epipelagic","Upper mesopelagic","Lower mesopelagic","Bathypelagic" ))SES_combined$indice <-factor( SES_combined$indice,levels =c("FRic","FDis","FDiv","FEve"),labels =c("Functional richness","Functional dispersion","Functional divergence","Functional evenness" ))ggplot(SES_combined, aes(x = depth_layer, y = SES, color = depth_layer)) +geom_point(size =3) +facet_wrap(~indice) +geom_hline(yintercept =c(1.96, -1.96), linetype ="dashed", color ="gray40", size=0.8) +scale_color_manual(values =c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +labs(x ="", y ="Standard Effect Size (SES)") +theme_light() +theme(axis.text.x =element_blank(),strip.text.x =element_text(size =14, color ="black"),strip.background =element_rect(fill ="white"),axis.title =element_text(size =13),axis.text =element_text(size =13)) +guides(col="none", fill="none")